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ABSTRACT 

We present here findings for C 18 depletion in eight starless cores in Taurus: TMC-2, L1498, 
L1512, L1489, L1517B, L1521E, L1495A-S, and L1544. We compare observations of the C 18 J=2-l 
transition taken with the ALMA prototype receiver on the Heinrich Hertz Submillimeter Telescope 
to results of radiative transfer modeling using RATRAN. We use temperature and density profiles 
calculated from dust continuum radiative transfer models to model the C 18 emission. We present 
modeling of three cores, TMC-2, L1489, and L1495A-S, which have not been modeled before and 
compare our results for the five cores with published models. We find that all of the cores but one, 
L1521E, are substantially depleted. We also find that varying the temperature profiles of these model 
cores has a discernable effect, and varying the central density has an even larger effect. Wc find no 
trends with depletion radius or depletion fraction with the density or temperature of these cores, 
suggesting that the physical structure alone is insufficient to fully constrain evolutionary state. We 
are able to place tighter constraints on the radius at which C 18 is depleted than the absolute fraction 
of depletion. As the timeline of chemical depletion depends sensitively on the fraction of depletion, 
this difficulty in constraining depletion fraction makes comparison with other timescales, such as the 
free-fall timescale, very difficult. 
Subject headings: s 



1. INTRODUCTION 

Understanding the transformation of a starless core 
into a protostellar object is the first step in under- 
standing the overall process of low-mass star formation. 
A key component in understanding low-mass star for- 
mation is determining the evolutionary states of star- 
less cores. The competing models for star formation 
- ambipolar diffusion and gravo-turbulent fragmenta- 
tion - ea c h hay e very different timescales (see , e.g. , 
Shu et all Il987t Klessen fc Ballesteros-ParedesI 120041: 
Mac Low fc Klessenl 120041: lEvanslll999t lHartmann et al.l 
20011) . 

Ambipolar diffusion posits that collapse of a starless 
core is slowed by neutrals slip ping relative to i ons and 
the decay of magnetic fields (|Shu et al.l 119871 ). The 
ambipolar diffusion timescale is approximately 5 x 10 6 
years, assuming an ionization fraction of 10~ 7 . Gravo- 
turbulent fragmentation proposes that large-scale tur- 
bulence causes fragmentation of individual cores, which 
then collapse on a free-fall timescale. For H2 densities of 
5 x 10 5 cm 3 , that implies a timescale of 5 x 10 4 years. 

One approach is to use chemical evolution, or the pro- 
cess by which the chemistry of a core changes over time, 
as another piece of information to distinguish between 
these theories. In the cold (< 10K), well-shielded inte- 
riors of dense starless cores, significant freeze-out (de- 
pletion) of molecular species s uch as CO, C S, CC S, 
and HCO+ are observed (e.g., iTafalla et all (|2004aD L 
The more evolved cores are expected to have more 
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freeze- out, or deplet ion (Tafall a~et al.l 1200 4a). For ex- 
ample, iCrapsi et all ((2007D found that in L1544, 93% 
of the C 17 has frozen out in the center of this very 
centrally condensed core (central density of 8 x 10 5 
cm -3 ). Another example is the work of Redman et al. 
(12002ft toward L1689B (central density of 2.8 x 10 5 cm" 3 
(jBacmann et~aIll2000D ) where 90% of CO is depleted af- 
ter 4 x 10 4 years, or half the free-fall timescale at that 
density. However, this is a lower limit, as they assume a 
sticking coefficient of unity, and a constant density over 
time (which, as they point out, is unreasonable - the den- 
sity would increase over time). 

The state-of-the-art theoretical studies are chemody- 
namical models that use a chemical network at each 
time step during the collapse of a starless core s to calcu- 
late a bundances. For instance, the models of iLee et al.1 
(2004) shows significant chemical depletion happening at 
roughly 10 6 years. It is therefore theoretically possible 
that the depletion profile can be used to distinguish be- 
tween ambipolar diffusion and gravo-turbulent fragmen- 
tation. 

In this paper, we map the J=2-l transition of C 18 to- 
ward a sample of eight starless cores in the nearby Taurus 
molecular cloud. By comparing observed integrated in- 
tensities with intensities obtained from radiative transfer 
modeling, we place constraints on the amount of deple- 
tion in terms of how much CO is depleted (fraction of de- 
pletion fd), and where this occurs (the depletion radius 
rd )■ Our work is unique because we use a larger sam- 
ple set (8 cores) , all located in the same molecular cloud 
and are then able to compare the amount of depletion in 
these cores to evolutionary parameters, such as central 
density. We do not assume a constant T g =10 K profile 
for our models, as has been done in previous work, and 
investigate degeneracies in choice of temperature versus 
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choice in abundance. We also investigate possible effects 
of the gas temperature and the dust temperature not 
being equal. Unless otherwise specified by Tdust> tem- 
peratures referenced are the kinetic gas temperature, T g 

This paper is organized as follows. In § 2 we 
discuss our sample selection and motivation. We 
also describe our observations and modeling method- 
ology, including the implementatio n of the RATRAN 
(|Hogerheijde fc van der Takl (|2000f ) ) radiative transfer 
code. In § 3 we test our modeling procedure, including 
potential degeneracies of various parameters, and present 
results for each individual core. In § 4 we present our 
analysis and in § 5 we give our conclusions. 

2. MOTIVATION AND PROCEDURES 
2.1. Sample Selection and Motivation 

Our sample consists of cores in Taurus, at approxi- 
mately 140 pc (|Loinard et alj|2005l ). so they are easy to 
compare at the same spatial scale. We chose this partic- 
ular set because they span over two orders of magnitude 
in central density (see Table 1). Figure Q] shows obser- 
vations of the dust continuum emission at 850 /jm taken 
with SCUBA, as well as C 18 2-1 integrated intensity 
maps. These cores have at least some rough central dust 
continuum concentration. 

Given a constant temperature profile (see § 3.2 for 
how this changes if temperature varies), and a constant 
[C 18 0]/[H 2 ], one would expect to observe more intense 
C 18 emission in the center that at outer radii. For 
many of these cores, that is not what is observed, as 
seen in Figure [T] Given the high central densities of these 
cores, the C 18 could be depleted by freeze-out onto dust 
grains. The greater the density, the more freeze-out one 
would expect, so this freeze-out would occur pre feren- 
tially in the centers of the cores (|Goldsmithl 12001 . We 
shall model this effect with radiative transfer. 

2.2. Observations 

We took observations of the J=2-l transition of C 18 
(219.56036 GHz) during February 5 - 8, 2007, using 
the sideband-separating ALMA 1mm (Band 6) proto- 
type receiver on the 10-m Heinrich Hertz Submillimetcr 
Telescope (HHT). The beam size, (full- width half max; 
FWHM) at the observed frequencies is 34". The main 
beam has been measured to be Gaussian to detectable 
levels using maps of Mars (J. Bieging, private commu- 
nication, 2008). For the ALMA prototype receiver, the 
sidelobes are insignificant, as those are down by at least 
20dB. The C 18 line was tuned to the lower sideband 
with a 6.0 GHz IF. Receiver noise temperatures were be- 
tween 90 — 100 K with system temperatures < 250 K. 
Typical sideband rejections were greater than 13 dB and 
were ignored in the subsequent calibration. The ALMA 
1mm prototype receiver is both sideband separating and 
dual polarization resulting in 4 intermediate frequency 
outputs (Vpol LSB, Hpol LSB, Vpol USB, and Hpol 
USB) to the spectrometers. Unfortunately, we were only 
able to use one polarization (Vpol) of the lower side- 
band since there was only one high spectral resolution 
backend {Av < 100 kHz) available in February 2007. A 
Chirp Transform Spectrometer with 32 kHz spectral res- 
olution was used for all C 18 observations. The HHT 



uses the standard Chop per Wheel Calibration Method 
(|Penzias fc Burrushl973f h The data were baselined and 
calibrated using the CLASS reduction package. Mapping 
observations were obtained on a 10" grid covering 2' x 2'. 
Each position was observed for 2 minutes of integration 
time in position-switching mode. OFF positions were 
checked for C 18 emission. Calibration scans of Saturn 
were made every two hours after pointing and focusing 
the telescope. The final reduced spectra were placed on 
the T m b = T^/r] m b scale with an average main beam ef- 
ficiency of r} m b = 0.74. 

Rms errors for the observed integrated intensity points 
generally ranged from 0.1 to 0.2, the rms error of the 
observed position with the largest 07 for each core, is 
given in Table 1. Errors are defined as: 

CTj = y/AVli ne AV c nannel 

* Ot (1) 

where AVu ne is the velocity interval over the whole 
line, from line wing to line wing. AV c h anne i is the channel 
width, of 0.157 km/sec, and <jt is the baseline rms error, 
assuming a linear baseline. 

2.3. Modeling 

A radiative transfer calculation is appropriate for this 
work, as computing a column density from observations 
requires the assumption of a single excitation temper- 
ature throughout the core, which is incorrect. Also, 
because we are in non-LTE, the kinetic temperature is 
not equal to the excitation temperature along the entire 
line of sight. Therefore, we use the 1-D version of the 
Monte Carlo radiative transfer code , RAT RAN, devel- 
oped by Hogcrhei ide &; van der Takl (2000), to generate 
radial profiles of C 18 emission. RATRAN uses Accel- 
erated Lambda Iteration to efficiently reach convergence 
between the level populations and calculations of J v , the 
mean intensity of the radiation field. 

RATRAN inputs for our models include: 

1. C 18 molecular collisi on rate file from the 
LAMBDA database (see ISchoier et all (l2005l )L 

2. Dust grain properties: We use OH5, see 
lOssenkopf fc Henninel (|1994fl . Table 1, column 5. 
These opacities typically give the best fit in star- 
less core dust continuum radiative transfer models 
(see Shirley et al. 2005). 

3. Radial points: We used 100 radial shells, from 
20 AU to 30,000 AU, with equal logarithmic spac- 
ing. 

4. H2 density profiles: The density profiles are con- 
strained by the Shirley & O'Malia du st continuum 
radiat ive transfer models (see also iShirlev et al.1 
(2005) for methodology). The models assume 
isothermal, static Bon nor-Ebert spheres with var- 
ious central densities (Bo nno r 119561: lEbertl 119551) 
which are derived under assumptions of pressure- 
bounded hydrostatic equilibrium. See Figure [2] for 
our density profiles. While the B-E spheres are 
isothermal, their resulting temperature profile from 
the dust continuum radia tive transfer calcula tion is 
non-isothermal; however. lEvans et al.l ()200lD found 
that the correction to n(r) for the non-isothermal 



3 



case is neglig i ble fo r low-mass starless cores. See 
IShirlev et alT(t2005h for further discussion. 

Our central densities are not identical to those used 
i n oth er published work, specifically in I Young et al.l 
(2004), because Shirley & O'Malia use a finer cen- 
tral density grid in the modeling. W e also assume 
a H2 ortho-para ratio of 3:1, as in iFlower et al.l 
(2006), and find that our results do not depend 
strongly on o:p ratio. 

Running these density profiles through RATRAN 
requires choosing a cutoff radius, beyond which we 
assume H2 (and, by extension C 18 0) has its stan- 
dard ISM value (see below). The cutoff radius is 
necessary for modeling; it is not physically moti- 
vated. In reality, the edge of a core is difficult 
to define, and there is more likely to be a grad- 
ual decline in density versus a sharp discontinu- 
ity between the core and surrounding ISM, we in- 
vestigated how sensitive our results were to this 
choice. This was esp e cially important given that 
work by lOnishi et al.l (|1996l) found large amounts 
of C 18 in the ambient Taurus molecular cloud, 
not necessarily associated with starless cores. We 
chose 30,000 AU as a reasonable outer radius for 
our models, co nsistent with the .1 pc typical core 
size given by iBergin fc Tafallal (120071 ) . For the 
C ls O observations, we mapped only to ±60" (cor- 
responding to 8400 AU at 140 pc) from the dust 
continuum peak, as the dust maps were restricted 
by chopping (9 c hop/2 = 60"). Figure |3] shows that 
most of the contribution to the total column den- 
sity occurs in the center parts of the cloud, sug- 
gesting our choice of outer radius is reasonable. 

5. C ls O abundance profiles: We assume a step func- 
tion, wherein C 18 is depleted by some fraction fd 
at some radius r^- We define: 

h = ([C 18 Q]/[H 2 ]) JSM 

([C 18 0]/[H2])mterior to r d 

We vary fd from 1.0 (no depletion) to 10 3 in steps 
of 2 n , with n varying from to 10. We vary rd 
from 1,250 AU to 30,000 AU in steps of 1,250 AU. 

Exterior to r d , [C 18 0]/[H 2 ] has its standard ISM 
abundance, which we define in this work as 1.6 x 
10~ 7 , hcraftcr referred to as the "canonical abun- 
dance" . Thi s is sl ightly lower than work by 
iFrerking etH] (|1982f ). their 1.7 x 10" 7 vs our 
1.6 x 10~ 7 ; we adjust it t o be consistent with work 
by I Wilson fc Roodl ()1994f ). who find a slightly lower 
18 0/ 17 ratio. Using Xr fitting, we find an rd and 
fd that best match both the observed spectra for 
the central region and the observed integrated in- 
tensity of C 18 emission across the core. Previ- 
ou s work on two of our cores, L1498 and L1517B, 
bv lTafalla et al.l (|2004a|) has shown that the abun- 
dance drop is quite sharp and t herefore th a t a st ep 
function is reasonable (see also iLee et al.l ((2003)) . 
See Figure |4] for our C 18 density profiles. 

6. Temperature profile of C 18 0: The temperature pro- 
files are self-consistently calculated from dust con- 



tinuum radiative transfer models of the cores. The 
dust continuum emission from the eight cores in 
this study were modeled as part of a larger sample 
of starelss cores that were observed with SCUBA 
(Shirley & O'Malia, in prep). A one dimensional 
dust continuum ra diative transfer code, CSDUST3, 
lEgan et"aT1 (|1988l) is used to match submillime- 
ter intensity profiles and spectral energy distribu- 
tions (SEDs). The heating in the dust continuum 
models is supplied by the Interstellar Radiation 
Field (ISRF) which is allowed to vary in overall 
strength and the degree of extinction (see Shirley & 
O'Malia, in prep.) For our purposes in RATRAN, 
we approximate the T(r) profile as the same for 
the gas kinetic t emperature and th e dust temper- 
ature. Work bv lYoung et al.1 (|2004l ) has indicated 
that this assumption holds true for densities above 
10 4 in L1498, above 10 3 in L1512, and above 5 x 10 4 
in L1544. Discussion of the implications of these 
assumptions can be found in §3.2. Figure [5] shows 
the temperature profile used in each core. 

7. Velocity profile: We use a static velocity field with 
turbulent and thermal line broadening described by 
the Doppler b parameter = 0.2, since that quantity 
fits the observed line widths. Late stages of collapse 
would produce li ne widths bro ader than those ob- 
served. Work by I M vers! (|2005f ) has indicated that 
the shape of the B-E sphere density profiles do not 
change substantially even in the collapsing case, 
until late in the collapse history. 

After the model has been run thro ugh RATRAN, we 
use the Miriad image analysis package (|Sault et al.lll995l) 
to convolve the outputs with our 34" beam size, and inte- 
grate over velocity to determine the integrated intensity 
in K km/s. 

Before we began any of our modeling, we first at- 
tempt to replicate previously published results, as a 
check on our procedures for running the radiative trans- 
fer code. We were graciously provided density and tem- 
perature profiles by Mario Tafalla from work on LI 52 IE 
IjTafalla fc Santiago! l200l . See Fig [6] for those profiles, 
as compared to our inputs for the same core. We were 
able to replicate their published best fit curves through 
our analysis, even though t hey used a different radiative 
transfer code (jBernesll 19791 ) . This indicates that any dif- 
ferences in our results are based on differences in real 
physical inputs, such as density and temperature, and 
not differences in modeling procedures. L1521E results 
are further discussed in § 3.3. 

3. DISCUSSION 

We begin our modeling by using model temperature 
and density profiles as discussed above, and then varying 
rd and fd, as in previous work to find best fits, rd and fd 
are considered free parameters and are fully explored as 
such. As these parameters represent the amount of deple- 
tion in each core, which are then used to make chemical 
evolution time scale arguments, it is important to under- 
stand how sensitive the results are to the assumptions 
that go into their calculation. We explore how much 
these calculated values of rd and fd would change if the 
RATRAN inputs changed. Of the seven RATRAN in- 
puts described in § 2.3, we explore factors influencing the 
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canonical abundance and temperature profiles of C 18 0. 
Inputs for molecular collision rates, dust grain proper- 
ties, and velocity profiles are already as well constrained 
as they are likely to be for now, and there is good agree- 
ment in the literature, so we do not explore those further. 
H2 density profiles are well-constrained by dust contin- 
uum modeling, and choice of cutoff radius is of little con- 
sequence to the model column densities. 

We choose to explore the canonical abundance of C 18 
because there is a large amount of disagreement in the 
literature over this value, and as we ourselves have no 
good way of independently verifying its correct value. 
We choose also to explore the sensitivity to the temper- 
ature profile, because while we self-consistently calculate 
a varying profile, other work assumes a constant profile. 
Furthermore, there are physically motivated reasons to 
believe that choosing different values of canonical abun- 
dance and/or temperature than described in the previ- 
ous section would yield different best-fit ra and fd- The 
point of this analysis is not to find an absolute best fit 
for this multi-degree parameter space, since input data is 
too limited to make such a finding meaningful. The point 
is to see how much caution is needed - if we use and 
fd as proxies for chemical evolution, and if these proxies 
are extremely sensitive to our input assumptions, then 
our work and future work needs to be mindful of this. 

3.1. Variations in the ISM C 18 O abundance 

As discussed in § 2.3, we assume that, in the absence of 
freeze-out, the ratio of [C 18 0]/[H 2 ] is fixed for all of our 
cores at the canonical abundance ofl.6xl0~ 7 . However, 
not everyone agrees that this is the correct value to use. 
A comparison of published work in the field shows that 
values of the canonical [C 18 0]/[H2] a bundance differs by 
a factor of seven, fr om 0.7 x 10~ 7 (Tafalla fc Santiago! 
12001 to 1.5 x 1 0~ 7 dTafalla et alJl2004ah to 4.8 x 10" 7 
( Lee et aLll2003|) . It may be that more than one of these 
values is correct, as the canonical abundance could even 
vary across the Taurus region. This may be a reasonable 
possibility in Taurus, given the large scale structure of 
filaments, cavities, and rings and subsequent variatio ns 
in the radiation field found by [Goldsmit h et al.l ((2008). 

Because of the wide range of published values for this 
particular input parameter, we choose to explore it in 
further detail by varying the canonical abundance and 
seeing how that affected fd and r<j. A full exploration of 
all possible values of canonical abundance is beyond the 
scope of this work, but we felt it necessary to at least 
gain a quick understanding of the nature of the impact. 
We ran models with all of the same input parameters de- 
scribed in § 2.3, but for half the canonical abundance, as 
well as twice the canonical abundance. In changing the 
canonical abundance, some of the cores were best fit by 
entirely different values of fd and Td than for an abun- 
dance ratio of 1.6 x 10 -7 (see next section for further 
discussion). To our knowledge even this cursory explo- 
ration of variances in canonical abundance has not been 
done before. Further constraints on canonical abundance 
are needed to be confident that fd and values generat- 
ing lowest Xr for choice of one canonical abundance have 
physical meaning. More details can be found in § 3.3. 

3.2. Sensitivity to Temperature Profiles 



We also explored how sensitive best fit fd and val- 
ues were to choice of temperature values. We started this 
investigation by using all of the same input parameters 
as before (see § 2), but changing the temperature pro- 
files. We explored both temperature profiles that vary 
and those remaining constant across the core. Our ob- 
jective was not, in this way, to constrain the temperature 
profiles of these cores but to understand how sensitive re- 
sults are to this parameter. We believe we are able to 
constrain well the temperature profiles from dust con- 
tinuum radiative transfer modeling, described in § 2.3. 
Additionally, a varying temperature profile, as used in 
our work, makes physical sense: the areas near the cen- 
ter of the core, being shielded from the interstellar ra- 
diation field, should b e colder than those on the outside 
(jZucconi et al.1 (|2000| )). Dust continuum modeling of our 
sample of cores suggests that while the temperature vari- 
ation is only a few Kelvin, this means that the outer ra- 
dius of the core is up to 30% warmer than the center 
(Figure [5]). 

However, much of the other work in this field has as- 
sumed that the temperature remains constant through- 
out the core, with a typical value of T=10 K for all radii. 
Given that our work does not make this assumption, we 
explore how temperature profiles could affect which fd 
and rd best fit the observational data. It is, of course, 
the combination of temperature and density that produce 
the observed intensity. If a core is cooler than assumed, 
then one would require a higher density to match the ob- 
servations, and vice versa. This would lead one to believe 
a core is less depleted than if one started with a higher 
temperature. 

We ran a representative set of models for B-E spheres 
with central densities of 3 x 10 5 and 2.5 x 10 4 (see Figures 
[7]and|S]). These Figures show the C 18 2-1 intensity ver- 
sus radius for various combinations of r<2 and fd, using 
various temperature profiles. We plotted three constant 
temperature profiles: T—7 K, T=10 K, and T=13 K, as 
well as a temperature profile that varied with radius be- 
tween 8 and 13 K. For Fig [TJ note how close the red line 
(temp=10 K) is to the black line (the varying tempera- 
ture profile). Since the varying temperature profile goes 
from 8-12 K, it seems one can simply take the average of 
those numbers and use that as the single temperature. 
Indeed, the shapes of all the cu rves are fairly similar . 
This is consistent with work by I Tafalla et al.l ((2004a) , 
who found that modeling with isothermal profiles can 
well reproduce the data. However, the exact choice of 
that isothermal temperature is important. From Figures 
[7] and 13 one can see that not all cores can be modeled 
at T=10 K. One would need further information, such as 
self-consistent temperature modeling from dust emission, 
as in this work, to know which isothermal temperature 
to use. Figure [5] shows the temperature profiles we adopt 
for each core in the rest of this work. 

Another common assumption is that the gas and dust 
temperatures are well coupled (we make that assump- 
tion in our work here). At high densities this is true, 
as collisions between the gas and dust temperatures re- 
sult in equal kinetic gas and dust temperatures. At low 
densities, however, as in those found on the outside of 
the core , this is not necessa rily the case. Indeed, recent 
work bv lYoung et al.l ()2004f ) has shown that gas and dust 



5 



temperatures can vary substantially in those conditions. 
However, precisely because Tj, ^ T gas only on the outside 
of the core, in the lowest density areas that contribute 
least to the column density (see earlier section), we do 
not expect this to be a large effect. We were graciously 
provided gas and dust temperature profiles by Young et 
al for three cores for which they did detailed energetics 
calculations: L1498, L1512, and 1544. We used our same 
model density profiles but their calculated gas and dust 
temperature. See § 3.3 for further discussion. 

We will see in section § 3.3 that in many cases we find 
different fd and r d values for many of the cores than 
previously published. Based on the above analysis, we 
believe some of this discrepancy is likely due to our tem- 
perature profiles being different from those used in other 
work. As another quick check, we simultaneously varied 
both the canonical abundance and the temperature for 
a small slice of parameter space. We found that there 
can also be a degeneracy in the abundance and the tem- 
perature used, as in Fig [9j which worsens the problem. 
One can well fit the data by assuming a higher canonical 
abundance and a lower temperature, or a lower canonical 
abundance with a higher temperature. 

3.3. Summary of Findings for Each Core 

Figure [10] plots the C 18 integrated intensity \r f° r 
depletion profiles with combinations of fd and rd- Due 
to the high scatter in the observed intensity points, we 
azimuthally averaged our observed points before compar- 
ing them to our 1-D model. We divided the data into 10" 
annuli around the center core, and averaged the observed 
and model points over that range. We calculated a Xr as 
follows: 



(J-obs Imorf) 



a bs * (number of annuli — 1) 



(3) 



We define l b s as the mean intensity of each of the 
observed points, for each 10" annuli. We define l mo d 
as the intensity of the middle model point in the same 
annuli. We define a oos as the variance in the intensity 
of the observed points for each annuli. For annuli with 
only one observed points (for some of the cores, there 
was only one observed data point in the 0-10" annuli), 
we define the variance as simply the error on that one 
point. Errors for these individual points are defined in 
equation 1. We have six annuli, going from to 60". 

For TMC-2, L1521E, L1495A-S, and L1544, no com- 
bination of fd and rd with a canonical abundance gave 
& Xr <!• However, as seen in Figure [TTl increasing the 
canonical abundance by a factor of two was able to give 
X r <1 values under 1, except for L1544 (discussed be- 
low). L1512 and L1517B had combinations of fd and rd 
that gave Xr <1 f° r both the canonical abundance and 
twice the canonical abundance, however using twice the 
canonical abundance ultimately gave lower Xr' s - None of 
the cores had lower Xr when the canonical abundance was 
halved, so we do not favor this possibility. In addition to 
these Xr maps, we also calculate Xr of the line profiles for 
the central 10" regions of these cores, as shown in Figure 
[T2l Figure [T3] shows observations of the central regions 
with best fit models (see Table 2). Both Xr calculations 
have significant overlap between the best fit values of 



fd and rd- While all of these models use the temper- 
ature profile calculated from dust continuum radiative 
transfer models, Figure [T4l plots the integrated intensity 
Xr for three cores using the temperature profil es from 
the energetics calculations of lYoung et al.l (|2004| ). where 
Td =/= T gas for all radii. This gives very similar results 
to Figure [TQ1 consistent with the conclusion that using 
uncoupled gas and dust temperatures does not produce 
a substantial effect on the C 18 J=2 — 1 emission. 

Inspection of Figures [TO] and QTJ shows that there is 
some degeneracy between exactly where the depletion 
occurs and by exactly how much. In many cases it is 
easier to constrain rd than fd- These figures also show 
the difficulty in constraining fd, as, & range of fd values 
all give Xr < 1 • Fig HH shows the observations for each 
core, as well as models with various Xr values. Despite 
this degeneracy, it is clear that L1521E has low depletion 
while cores such as L1512 and L1517B, for example, have 
higher depletion. The best fit fd and rd values for each 
source are listed in Table 2. Details for each core are 
given below. 

1. TMC2: this is a difficult core to model, since, based 
on dust emission (see Figure[TJ), it is not as centrally 
condensed, and is non-azimuthally symmetric. In 
the case of canonical abundance, it is best matched 
by a model profile with no depletion. However, 
even this profile has a fairly high Xr- I n order to 
get the lowest x% we had to double the C 18 ISM 
abundance. The lowest x r f° r TMC-2 is a depleted 
model with twice the canonical abundance. 

2. L1498: We find this core is well-modeled with 
canonical abundance of C ls O with significant de- 
pletion. We find s imilar depletion to work by 
I Willacv et al.l ()1998l ) , which sets a l ower limit on 
the dep letion fraction at 8. Work by iTafalla et al.l 
(|2004al ) found C 18 to be entirely depleted {f d > 
1000) inside a radius of 9,940 AU. We find a r d of 
7,500 AU but only a factor of 16 depletion. How- 
ever, as shown in FigfTU l we find that our findings 
and those of lTafalla et al. (2004a) to have very sim- 
ilar Xr- The slight discrepancy may be due to our 
azimut hally averaging our data, while ITafalla et al~l 
(|2004aft did not. 

3. L1512: The best match for this core is twice the 
canonical abundance of C 18 with sign ificant de- 
pletio n. This is consistent with work bv lLee et al.l 
(2003). They give a greater rd than this work 
(15,000 AU instead of this work's 10,000 AU) but a 
lower fd (25 versus this work's 256). As in the case 
of L1498, Fig [TTl shows both findings to have simi- 
lar Xr ' s i an d could be exp lained by our azi muthally 
averaging our data, while iLee etaf! (|2003l ) did not. 

4. L1489: This core is well-modeled with canonical 
abundance of C 18 with significant depletion. 

5. L1517B: Our results for L1 517B show s i gnifican t 
depletion, as does work by ITafalla et al.l ( 2004bfl . 
Their work, however, found a best fit (to their ob- 
servations) with f d > 1000 and r d = 11,620 AU, 
which is more depletion than our work (see Table 
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2), perhaps because they use a T=10 K tempera- 
ture profile, higher than ours. 

6. L1521E: Our best fit for this core is twice 
our canonical abundance with 7^=2,500 AU and 
/d=64. This represent s little or no depleti o n, con - 
sistent with work by iTafalla fc Santiago! (|2004T ). 
and suggests an ear ly evolutiona r y stat e, consistent 
with work done by lHirota et al.l (|2002l ). Our exact 
findings for and fd , are not identical to those 
in this work, however, as lTafalla fc Santiago! (j2004) 
finds L1521E best fit by no de pletion at all. How- 
ever, ITafalla fc Santiago! (|200l use a different den- 
sity profile, and a different canonical abundance of 
[C 18 0]/[H 2 ] (ours is 1.6 x 1CT 7 , versus 0.7 x 1CT 7 ), 
and assumes a constant temper ature T=10 K at all 
radii (Figure U3). Additionally. ITafalla fc Santiago! 
()2004l ) observational values for the C 18 2-1 inte- 
grated intensity are 30% higher than ours (Figure 
[HI). Discussions with Tafalla have indicated that 
some of this discrepancy may be the result of dif- 
ferent beam sizes, as the FWHM of his observations 
is 12" versus our 34", and systematic calibration er- 
rors (the IRAM observations were double sideband, 
while the HHT observations were single sideband). 

7. L1495A-S: Very high Xr, n °t well matched by any 
Td and fd with canonical abundance. Similar to 
TMC-2 in that it is best matched by a profile with 
no depletion if one assumes a canonical abundance. 
The lowest Xr value, however, comes from doubling 
the abundance and assuming significant depletion. 

8. L1544: Our results show this core to be best fit 
by using twice the canonical abundance but still 
significantly depleted, which is consistent with the 
conclusion of Lee et al (2004). They find an fd of 25 
at an of 9,000 AU. This is more depletion than 
our findings, fd of 16 at an Td of 7,500 AU, however 
they a re using a higher ca nonical abundance, 4.8 x 
10 -7 . iCaselli et al.1 (|1999| ) found this core to have 
an abundance drop by a factor of approximately 
10 at a radius of 6500 AU. Our Fig ITT1 shows this 
result to have a higher \r than our best fits, but 
still fa irly consistent with our results. lAikawa et al.l 
(pOOlh find a central hole of 4,000 AU, inside of 
which no CO is present. This is more depletion 
than our findings, but their work assumes a higher 
H2 number density than ours, 1.5 x 10 6 cm -3 . Our 
best fit value still has a Xr greater than one, as seen 
in Table 2. This is likely due to L1544 being not 
round, but inspection of Figure [15] shows that we 
were able to find a combination of fd and that, 
at least by eye, appears to be a good fit. 

4. ANALYSIS 

We find that all cores but one, LI 52 IE, are best fit by 
models that include significant depletion. Our conclusion 
that all but one of these cores is significantly depleted 
holds even taking into account reasonable variances in 
outer radius, temperature, and canonical abundance. In 
no case did we find varying any of these parameters sub- 
stantially affected the results. This work is unique in 
that we have tested a more complete parameter space of 



r d, fd, and C 18 canonical abundance than previously 
published work. While outer radius does not impact 
our results, choices of temperature profile and canoni- 
cal abundance can. Even fixing these parameters, we 
cannot, however, robustly distinguish between a factor 
of 10 depiction and a factor of 1000 in fd, as the \r are 
very similar (e.g., . Figures Q~5l ITUl and ITTj) . In con- 
trast, the depletion radius, rd, for some of our cores is 
more strongly constrained than fd] however, values 
for many of these cores are separated by little more than 
a beam width, a nd four of t he cor es have the same value 
for rd- Work bv lLee e~ al. (2003() also found fd hard to 
constrain. 

At the distance of Taurus, information in the HHT 
maps contain information on scales larger than m b/2 — 
2100 AU. Greater spatial resolution would better con- 
strain rd- Improved resolution may even help constrain 
fd- For instance, future observations with the Large 
Millimeter Telescope (LMT) 50-m operating at 1.3mm 
would increase our resolution by a factor of 5. The 
difficulty in determining the absolute value of fd lim- 
its our ability calculate an accurate chemical depletion 
time scale. The rate of freeze-out, or the chemical deple- 
tion timescale of CO, is given in iRawlings et al.l (119921) 
and | Redman et ail (|2002f ) . The version in Redma n et al.l 
gives h C o oc n C o- The solution of this differen- 
tial equation implies that the amount of gas-phase CO 
decreases exponentially with time for a constant grain 
temperature, nco The inability to tell the dif- 

ference between /d=10 vs /d=1000 depletion translates 
into a factor of 4 diffe r ence i n the timescales. As pointed 
out in lRedman et al.1 |2002), for a core with H 2 density 
of 2.8 x 10 5 (near the median central density of cores 
modeled by Shirley & O'Malia, in prep.), 90% of the 
CO is depleted at 43,000 years, or approximately half 
the central density free-fall time scale. But if the CO is 
actually depleted by 99%, this brings it to the free-fall 
time scale, and 99.9% depletion would be longer than 
the free-fall time scale. The difference between /d=10 vs 
/d=1000 could mean the difference between a core be- 
ing much older than the free-fall time scale, hence dis- 
favoring gravo-turbulent fragmentation as the dominant 
core evolution scenario, or being younger than the free- 
fall time scale, and consistent with gravo-turbulent col- 
lapse timescales. 

Better constraints on the amount of depletion 
will be important for constraining ambipolar diffu- 
sion timescales, which d epend on ionization fraction. 
ICaselli fc Walmslevl (pOOlh found that freeze-out of C 18 
can be accompanied by depletion of metal ions, signifi- 
cantly lowering the ionization fraction. In the case of 
L1544, they find that ionization fraction can be 10~ 9 , 
which would put the ambipolar diffusion timescale on the 
order of the free-fall timescale. They argue that better 
constraints on amount of CO depletion could place con- 
straints on the ionization fraction, so these issues could 
be coupled. This timescale could be shorter than the 
chemical evolution time scale (depending on the amount 
of depletion, as discussed earlier), which would imply 
ambipolar diffusion is not the dominant process. 

In Figure [TU we compare the best Xr model parame- 
ters to the physical parameters of the cores (central den- 
sity, central temperature, as well as density and temper- 
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ature at best fit values). The plotted points indicate 
the best fit fd and (minimum x 2 .), as given in Table 
2. Error bars show the range of fd and Td that corre- 
spond to Xr < min(Xr)+l- Figure [151 shows that it is 
much more difficult to constrain fd than r^. We find no 
strong correlations for any of these parameters with r^, 
and none for fd, although large error bars on fd make 
teasing out trends difficult. One would expect the chem- 
ical evolution to track the dynamical evolution; as a core 
contracts, more CO would freeze out. Our work suggests 
something more complicated, as the cores with greater 
central density are no more likely to be significantly de- 
pleted than those with low density. Warmer cores are 
also no less likely to be depleted. For instance, L152fE 
has one of the lowest dust temperature profiles yet has 
the least depletion. 

Our findings are consistent with the conclusion that 
these cores may be evolving at different dynamical rates. 
There are several possible explanations for why cores in 
Taurus could be evolving at different rates. For exam- 
ple, that the magnetic field strengths are different across 
Taurus. Since the magnetic pressure oc B 2 , small changes 
in magnetic field would lead to large changes in pressure 
support. This would lengthen the collapse timescale until 
the cores become supercritical. Unfortunately, the mag- 
netic field strengths in Taurus starless cores are poorly 
constrained and the observed variation from dust polar- 
ization measurements is over an order of magnitude (e.g. , 
» IOmG inL1498 to > 100/^G in L1544 (|Kirk et aLll200a 
ICrutcher et aDl200l . 

Rates of evolution may also be affected by turbu- 
lence. For instance, the external pressure is an impor- 
tant fac tor in de t ermin ing the gravitational stability of 
the core iBonnorl (|1956l) . While the starless cores them- 
selves have very narrow linewidths that are not domi- 
nated by turbulence, the turbulence outside the cores 
may play an important role in the bounding pressure. 
The transition from thermal to tu rbulent medium ca n 
be quite sudden in dense cores ([Pineda et all I20T0T ). 
The amount of turbulent pressure varies across Tau- 
rus ((Goldsmit h et al.l [2008). This external pressure is 
in addition to the weight of clo ud material surround- 
ing the cores (jLada et al.1 [2008) . Observations of the 
surrounding turbulent environment (e.g. 13 CO maps, 
Goldsmit h et al.l ((2008)) compared with observations of 
the transition between the coherent core and more turbu- 
lent surrounding medium (e.g. NH3 maps with the new 
7-beam receiver at the Green Bank Telescope) are needed 
to elucidate the importance of the bounding pressure and 
the role of turbulence on core evolution. Further infor- 
mation on the effects of magnetic fields and turbu lence 
specific to Taurus can be found in lLee et al.1 fj2009f ) and 
iGoldsmith elaH (poM ). 

It is also possible that freeze-out is not as simple as we 
may expect, a s deso rption processes may be important. 
iRoberts et al.l (|2007f) investigated the effects of desorp- 
tion mechanisms such as X-rays heating, direct impact of 
cosmic rays, UV irradiation, and exothermic reactions on 
dust grain surfaces (e.g., formation of H2). They dismiss 
X-ray effects, as there are not many X-rays in these types 
of molecular clouds, but find that the other three mech- 
anisms can potentially be significant. Calculating the 
magnitude of these effects, however, depends on poorly 



constrained parameters, such as the composition of ice 
on the mantle and the size and sh ape of dust gra ins 

Worsening the problem is that IRoberts et al.l ([20071 ) 
are forced to use observations of CO freeze-out to con- 
strain these parameters. As a result, if there is a sub- 
stantial uncertainty between the amount of depletion 
as found in our work, then these parameters are even 
more poorly constrained. They conclude that freeze-out 
rates are still high (98% for a quiescent core of density 
10 5 ) even with these effe cts; however, they ad mit a large 
amount of uncertainty. lAikawa et al.l (|2003l ) have also 
done work on surface grain interactions, which could also 
help in understanding desorption processes. More work 
is needed in this area. 

Velocity structures and dynamical evolution may be 
tied to chemi cal evolut i on in more indirect ways, as 
pointed out in iLee et afl (|2005l ) . They compare chemical 
evolution in static and collapsing density profiles, and 
find that in prestellar cores, both give similar results. 
However, in a dynamic case, the material being pulled 
into the center of the core may be coming from an outer 
radius that has less depletion. A better understanding of 
these cores' dyn amics and env i ronm e nts would be ver y 
useful. Workbv lATkawa et al.l (|2005l ). ILee et al.l (|200l . 
and iLed ([20051 ) has made much progress on this topic. 
However, as pointed out by M. Tafalla (private commu- 
nication, May 9, 2010), the linewidths observed in many 
of these cores are inconsistent with arguments that cores 
such as L1521E are simply collapsing too quickly to have 
time for chemical freeze-out. 

5. CONCLUSIONS 

We find that all of our cores but one, L1521E, are best 
fit by models showing significant depletion, consistent 
with previous work in the field suggesting that unevolved 
cores are rare. This would imply that the depletion pro- 
cess happens very quickly. The lack of unevolved cores 
is likely an observational bias, as chemically unevolved 
cores may also be lower density than those in our sample 
set. Observations with Herschel and SCUBA2 may find 
many more lower density cores that may also be chemi- 
cally unevolved. 

We find that assuming an isothermal temperature pro- 
file can reproduce the observations, although one needs 
further information to know the appropriate temperature 
to use. This is likely due to the small temperature vari- 
ations seen in dust continuum radiative transfer models. 
T=10 K, however, rarely provides a good best match. 
We find assuming Td — T gas gives similar results to ener- 
getics calculations with Td ^ T gas . Our work also finds 
no correlation between the amount of depletion and cen- 
tral density or central temperature, nor do we find cor- 
relations between density or temperature at the radius 
of depletion. It is possible that these cores are evolving 
at different rates that may be affected by variations in 
the magnetic field strengths, turbulence, and bounding 
pressures across Taurus. 

Constraining the amount of depletion in these cores 
can be complicated by degeneracies between tempera- 
ture and canonical abundance and difficulty in placing 
strong constraints on fd- Constraining fd more strongly 
is necessary to compare chemical evolution time scales to 
free-fall time scales, an important step in distinguishing 
between the relative importance of ambipolar diffusion 
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and gravo-turbulent fragmentation in the evolution of 
starless cores. Ultimately, observations with better an- 
gular resolution, as well as additional understanding of 
the desorption processes, dynamics, and surrounding en- 
vironments of these cores, will help complete our under- 
standing of depletion and chemical evolution in starless 
cores. 
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TABLE 1 

Characteristics of These Cores 



Core 


RA (J2000) 


Dec 


V LSR 


Central density Tinner Touter 


a I(C ls O 2-1) 








km/sec 


cm 3 K K 


K*km/sec 



TMC-2 


04:32:45.5 


+24:25:07.7 


+6.30 


2 x 10 4 


10.7 


10.9 


0.180 


L1498 


04:10:53.0 


+25:10:08.6 


+7.78 


2 x 10 4 


10.2 


10.3 


0.194 


L1512 


05:04:08.6 


+32:43:24.6 


+7.10 


8 x 10 4 


8.2 


8.6 


0.157 


L1489 


04:04:47.6 


+26:19:17.9 


+6.85 


1 x 10 5 


11.7 


12.4 


0.166 


L1517B 


04:55:18.3 


+30:37:49.8 


+5.78 


2 x 10 5 


8.7 


9.4 


0.200 


L1521E 


04:29:14.9 


+26:13:56.6 


+6.90 


3 x 10 5 


7.6 


8.5 


0.187 


L1495A-S 


04:18:39.9 


+28:23:15.9 


+7.29 


4 x 10 5 


9.7 


12.8 


0.187 


L1544 


05:04:17.2 


+25:10:43.7 


+7.18 


8 x 10 5 


8.5 


11.4 


0.166 



TABLE 2 

Best Fit: rd, fd, and Abundances (as a multiple of the canonical abundance) 



Core 


ra 
AU 


fd 


Abundance 


Xr 


TMC-2 


5000 


512 


2 


0.094 


L1498 


7500 


16 


1 


0.063 


L1512 


10000 


256 


2 


0.13 


L1489 


7500 


1000 


1 


0.48 


L1517B 


10000 


32 


2 


0.087 


L1521E 


2500 


64 


2 


0.18 


L1495A-S 


7500 


4 


2 


0.42 


L1544 


7500 


16 


2 


1.5 



10 



850 fim C 18 J=2-l 850 fim C 18 3=2-1 




Fig. 1. — On left, 850 fim dust continuum emission, from SCUBA data, normalized (arbitrary units). On right, C O integrated intensity. 
Beam size indicated by solid circles. Contours are percentage of the peak for the 850 /imcontinuum, and value indicated for the C ls O 
integrated intensity. 
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Fig. 2. — H2 number density vs. radius profiles used for each of the cores. 
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Fig. 3. — We explore possible effects of choice of Router- We integrate the column density from r=0 (the center) to r=R ou t er (the outer 
radius, which we vary from zero to our model cutoff of 30,000 AU), then divide that value by the integrated column density from to 
30,000 AU, and plot that ratio as a function of R ou ter • We do this for three model cores with densities of 2 X 10 4 , 3 X 10 5 , and 8 X 10 s cm~ 3 , 
representing the full span of our sample's central density. 
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Fig. 4. — Density of C 18 vs. radius for representative combinations of (in AU) and fd- Overplotted (red dashed line) is the density 
of H2, just shifted down by the canonical abundance, [C ls O]/[H2]= 1.6 X 10~ 7 , so as to fit on the same plot. 
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Fig. 5. — Temperature vs. radius profiles from dust continuum radiative transfer models used for each of the cores. 
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Fig. 6. — For L1521E, comparison of temperature and density profiles used in Tafalla & Santiago (2004) to those in this work. The 
bottom panels compare observed and model values of integrated intensity (K km sec -1 ) of the C ls O J=2-l for both this work and theirs. 
In the bottom left panel, blue crosses are their observed points. The solid black line was made by using their temperature and density 
profiles, which we ran through RATRAN and our associated analysis. It appears to be identical to their best fit curve found in their Figure 
3. In the bottom right panel, the black points are our observed points. The red (upper) line is for no depletion, with our temperature and 
density profiles (twice the canonical abundance), run through RATRAN. The green (lower) line is our overall best fit with r<j=2,500 AU 
and /tj=64. For further discussion, see § 3.3. 






Fig. 7. — Intensity vs. radius for various combinations of and fd- This is for a theoretical B-E sphere with an example density of 
3 X 10 5 cm -3 . The bottom panel shows the temperature profiles used in the nine panes above, with matching linestyles. 
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Fig. 8. — Intensity vs. radius for various combinations of and f^. This is for a theoretical B-E sphere with a central density of 
2.5 X 10 4 cm -3 . The bottom panel shows the temperature profiles used in the nine panes above, with matching linestyles. 
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Fig. 9. — Illustration of modeling degeneracies present in L1521E, as an example. Note that our data points (black) can be fit well with 
either our temperature profile, which varies with radius and double the canonical abundance of C ls O, or with a constant temp profile of 
T=10 K and the regular abundance. Model points exactly match coordinates of observed points. 
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Fig. 10. — Xr maps of each of our 8 sources, assuming a canonical abundance, 1.6 X 10 -7 . Observed and model points have been 
azimuthally averaged in 10" annuli. The black regions represent Xr ^ ess than minimum(Xr)+l- For TMC2, L1521E, L1495AS, and L1544, 
no combination of and gives a Xr l ess than 1. Note that the a wide range of gives very similar Xr values, showing this value is 
poorly constrained. Blue diamonds represent the location of the lowest Xr> note they are in bottom left corner for TMC2 and L1521E. 
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Fig. 11. — Xr maps of each of our 8 sources, assuming twice the canonical abundance, or 3.2 X 10 -7 . Observed and model points have 
been azimuthally averaged in 10" annuli. TMC2, L1521E, and L1495AS, now have Xr l ess than 1. fd is still poorly constrained. Blue 
diamonds represent the location of the lowest Xr' note they are in the upper right corner for L1498 and L1489. 
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Fig. 12. — Xr °f the central regions, showing the models that best match the central spectrum, assuming canonical abundance. The 
contours (starting with the minimum Xr value and increasing by 1) are similar in shape to Figure [TOl but with higher Xr i as there is more 
noise in the line. L1495AS is not well fit due to the flat line shape, see Figure [13] Blue diamonds represent the location of the lowest Xri 
note they are in bottom left corner for TMC2 and L1521E. Note this figure is not azimuthally averaged since it is focused only on the 
central region of each core. 




Fig. 13. — Central spectra of each core, for both the observed values (black solid) and the model value (red dotted) with the lowest \ 
for the whole region. See Table 2. 




Fig. 14. — Xr maps of cores for which work by Young ct al. (2004) found that gas and dust temperatures are not the same. Data and 
model points have been azimuthally averaged in 10" bins; contour levels start with the minimum \r value. Blue diamonds represent the 
location of the lowest \r ■ 
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Fig. 15. — Observational data (in black), as well as the combination of and /<; that minimized the reduced \r ( m blue). See Table 
2. Errors in observed points are as defined in Equation 1. TMC2, L1521E, L1495A-S, L1517B, and L1544 are fitted with an abundance 
twice canonical. Plotted in black are the combination of and fd that give the highest \r> ^ ess than minimum(x?_)+l. Bla ck and blue 
lines are similar in most cases, illustrating the difficulty in constraining these quantities. See contour plots in Figs llOl and llll This shows 
the difficulty in constraining f^. For L1544, only one and fd combination gives \r < mm (Xr)+l- 
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Fig. 16. — We plot best fit and versus density and temperature at the centers of the cores. Error bars show the ranges of values 
that give Xr < rninimum(x^)+l. We find no strong correlation of or with these parameters, although error bars for are large. 
This shows that is harder to constrain than rj. The bottom two panes plot temperature and density for each core at the core's best fit 
values. These also show no strong correlation, as each core appears to be experiencing depletion at different temperatures and densities. 
L1544 does not appear in the figure, since no combination of and gave Xr <1 f° r that core. 



